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Abstract 

The Bhatnagar-Gross-Krook kinetic model of the Boltzmann equation is solved for 
the steady cylindrical Poiseuille flow fed by a constant gravity field. The solution is 
obtained as a perturbation expansion in powers of the field (through fourth order) 
and for a general class of repulsive potentials. The results, which are hardly sensitive 
to the interaction potential, suggest that the expansion is only asymptotic. A critical 
comparison with the profiles predicted by the Navier-Stokes equations shows that 
the latter fail over distances comparable to the mean free path. In particular, while 
the Navier-Stokes description predicts a monotonically decreasing temperature as 
one moves apart from the cylinder axis, the kinetic theory description shows that 
the temperature has a local minimum at the axis and reaches a maximum value 
at a distance of the order of the mean free path. Within that distance, the radial 
heat flows from the colder to the hotter points, in contrast to what is expected from 
the Fourier law. Furthermore, a longitudinal component of the heat flux exists in 
the absence of gradients along the longitudinal direction. Non-Newtonian effects, 
such as a non uniform hydrostatic pressure and normal stress differences, are also 
present. 
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1 Introduction 



The steady flow in a long channel or in a long tube of circular section un- 
der the action of a difference between the pressures imposed at the two ends, 
usually known as Poiseuille flow or Hagen-Poiseuille flow, is a typical text- 
book example in fluid dynamics [1] . In the last few years, a number of authors 
[2,4,3,5-11] have analyzed this problem with the channel geometry when the 
pressure difference is replaced by a constant external field g. Kadanoff et al. 
[2] have simulated this flow with the FHP lattice gas automaton [12] to con- 
firm the validity of a hydrodynamic description for lattice gas automata. For 
a dilute gas, Esposito et al. [4] have analyzed the solution of the Boltzmann 
equation in the Navier-Stokes limit. A generalized Navier-Stokes theory was 
seen to give a reasonable account of a fluid composed of molecules that possess 
spin when compared with molecular dynamics simulations [7]. Other studies 
[3,5,6,8-11], on the other hand, have focused on the breakdown of the contin- 
uum hydrodynamic predictions when the strength of the external field is not 
asymptotically small. In Ref. [3], an exact solution of the Bhatnagar-Gross- 
Krook (BGK) model kinetic equation was found for a particular value of the 
field strength. The general solution under the form of an expansion in pow- 
ers of g was considered in Ref. [5], where explicit expressions were derived 
to fifth order in the field. More recently, the general solution corresponding 
to the Boltzmann equation for Maxwell molecules has been derived to sec- 
ond order [8] and approximate solutions for hard spheres have been obtained 
from a Burnett description [10] and by means of moment methods [9,11]. 
The theoretical predictions of Ref. [5] have been confirmed at a qualitative 
and semi-quantitative level by numerical simulations of the Boltzmann equa- 
tion [6,10] and by molecular dynamics simulations [9]. The most surprising of 
those theoretical predictions is that the temperature profile exhibits a bimodal 
shape, namely a local minimum at the middle of the channel surrounded by 
two symmetric maxima at a distance of a few mean free paths. In contrast, the 
continuum hydrodynamic equations predict a temperature profile with a (flat) 
maximum at the middle. As a consequence, the Fourier law is dramatically vi- 
olated since in the slab enclosed by the two maxima the tranverse component 
of the heat flux is parallel (rather than anti-parallel) to the thermal gradient. 
Furthermore, non-Newtonian properties, such as normal stress differences and 
an effective shear viscosity depending on the hydrodynamic gradients, are also 
present. 

The goal of this paper is to carry out a kinetic theory description of the 
Poiseuille flow driven by an external force when the gas is inside a pipe. The 
reason is two-fold. First, the pipe geometry is much more realistic, and thus 
more worth studying, than the channel one. Second, it is important to test 
whether the failure of the continuum description to account for some of the 
qualitative features of the Poiseuille flow in a channel is not linked to that 
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particular geometry and extends to the pipe case as well. In this context, it 
is worthwhile noting that an exact solution of the Boltzmann equation (with 
Maxwell molecules) for the planar Fourier flow, which is valid for arbitrary 
values of the thermal gradient [13,14], does not extend to the cylindrical ge- 
ometry [13]. The results reported in this paper, on the other hand, confirm 
that the structure of the hydrodynamic and flux profiles in the pipe problem 
is quite similar to that of the channel problem. In particular, the temperature 
exhibits a non-monotonic behavior as one moves apart from the pipe axis. 
Nevertheless, the deviations from the continuum description are in general 
quantitatively smaller than in the channel case. 

The paper is organized as follows. The continuum hydrodynamic description is 
worked out in Sec. 2, both for the channel and cylindrical geometries. When the 
Navier-Stokes constitutive equations are inserted into the exact balance equa- 
tions for mass, momentum and energy, a closed set of coupled equations for 
the hydrodynamic fields (pressure, flow velocity and temperature) is obtained. 
The spatial dependence of the transport coefficients through the temperature 
is taken into account. Given the nonlinear character of the set of equations, its 
solution is expressed in powers of the external force. Section 3 is devoted to a 
summary of the results obtained from the kinetic theory description in the case 
of the channel Poiseuille flow and a critical comparison with the Navier-Stokes 
predictions is carried out. The original part of the paper is presented in Sec. 4, 
where the BGK kinetic equation is solved for the pipe Poiseuille problem by 
means of a perturbation expansion in powers of the force. Explicit expressions 
for the successive contributions to the velocity distribution function through 
fourth order in the force are derived. By velocity integration, the profiles of 
the hydrodynamic fields and their fluxes are then obtained. Since the content 
of Sec. 4 is rather technical, the discussion of the results is postponed to Sec. 5, 
where only terms through third order are considered. As in the planar case, the 
kinetic theory results strongly differ from the continuum theory expectations, 
especially in the case of the temperature profile. The breakdown of the Fourier 
law is characterized by an apparent thermal conductivity coefficient that only 
for large distances tends to the Navier-Stokes coefficient. Analogously, an ap- 
parent shear viscosity coefficient is introduced to monitor deviations from the 
Newton law. Finally, the main conclusions of the paper are presented in Sec. 
6. 



2 Hydrodynamic description 
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2. 1 Channel geometry 



Let us first consider the Poiseuille flow in the channel geometry. A fluid is 
enclosed between two infinite parallel plates normal to the y axis and located 
at y — ±i?, which are kept at rest. A constant external force per unit mass 
(e.g., gravity) g = gz is applied along a direction z parallel to the plates. 
We assume a laminar and incompressible regime, so in the steady state the 
physical quantities depend on the coordinate y only. The balance equations 
for momentum and energy become 

BP 
BP 

yz dy + dy~ U ' [6) 



where p is the mass density, u = u z z is the flow velocity, P is the pressure 
tensor and q is the heat flux. In the Newtonian description these fluxes are 
related to the hydrodynamic gradients by the Navier-Stokes (NS) constitutive 
equations. In this problem they read 

Pxx Pyy Pzz Pi (4) 

Pyz = _r? ^' (5) 
dT 

q z = o, (7) 



where p = |Tr P is the hydrostatic pressure, T is the temperature, and r\ and k 
are the shear viscosity and the thermal conductivity, respectively. Combining 
Eqs. (l)-(6), we get 

d du z 

d_ K dT = _ (duA 2 (1Q) 
dy dy \By J 
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Equation (9) gives a parabolic-like velocity profile, that is characteristic of the 
Poiseuille flow. The temperature profile has, according to Eq. (10), a quartic- 
like shape. Strictly speaking, these NS profiles are more complicated than just 
polynomials due to the temperature dependence of the transport coefficients. 
Since the hydrodynamic profiles must be symmetric with respect to the plane 
y = 0, their odd derivatives must vanish at y = 0. Thus, from Eqs. (9) and 
(10) we have 



d 2 u 7 



dy 2 



y=0 



d 2 T 



dy 2 



Po9 

Vo ' 



0, 



y=0 



d 4 T 



-2 



y=0 



P 2 o9 2 



(11) 

(12) 
(13) 



where the subscript denotes quantities evaluated at y — 0. According to 
Eqs. (11)— (13), the NS equations predict that the flow velocity as well as the 
temperature have a maximum at the middle layer y = 0. As we will see in 
the next Section, the kinetic theory description shows that the temperature 
actually exhibits a local minimum at y = 0, since d 2 T / dy 2 \ y=0 is a positive 
quantity (of order g 2 ). 

In order to get the hydrodynamic profiles from Eqs. (9) and (10), one needs to 
know the density and temperature dependence of the transport coefficients. 
To fix ideas, let us consider a dilute gas of Maxwell molecules, in which case 
p oc pT, rj oc T, k oc T [15]. Consequently, from (8)-(10) we get 

p(y) = po, (14) 
(y) = u - p ^f, (is) 

T(y) = T - -^^y\ (16) 
yy> 12^o 1 ^ 



u 



where y is an auxiliary space variable defined by dy = [T /T(y)]dy. From Eq. 
(16) the relationship between the true space variable y and the scaled quantity 
y can be found as 

y = y(l- p2 ° 9 \ f) . (17) 
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Equations (5) and (6) then imply that 

P yz (y) = pogy, (18) 

The solution of the fifth-degree equation (17), once inserted into Eqs. (15), 
(16), (18) and (19), yields the velocity, temperature and flux profiles predicted 
by the NS equations in the case of Maxwell molecules. To fourth order in the 
gravity field, the results are 

u , to) = ^_^( 1 + _^) +0(A (20 , 

^ r »-iC4 + ^ 4 H 96) < (21) 

P ^=^( 1 + ^f/) + °^ (22) 



It is interesting to note that y can be eliminated between Eqs. (15) and (16) 
to obtain the following nonequilibrium "equation of state": 

T = T -^(u -u z ) 2 , (24) 



which is independent of g. By equation of state we mean in this context a 
relationship holding locally among the hydrodynamic fields (p, u z and T) 
and that is independent of gravity (at least up to a certain order). Since the 
pressure is uniform in the NS description, it does not enter into Eq. (24). 
Interestingly enough, a quadratic dependence of the temperature on the flow 
velocity also appears in the case of the steady planar Couette flow [16] 



2.2 Cylindrical geometry 



Now we assume that the fluid is inside a straight tube of uniform circular 
section of radius R. Let the z axis be parallel to the pipe axis. As before, 
an external force per unit mass g = gz is applied to produce a flow field. In 
the (laminar) steady state all the physically relevant quantities depend only 
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on the distance r = (x 2 + y 2 ) 1 ^ 2 from the axis. In the case of this cylindrical 
geometry, the exact balance equations become 



d_ 

dr 



(vP rr ) P<f>4>i 



d 

T ~ X ~dr ^ Prz) = P9, 
Ou 



(25) 
(26) 
(27) 



These equations constitute the cylindrical counterpart of Eqs. (l)-(3). In gen- 
eral, the relationships between the cylindrical and Cartesian components of a 
vector A and a tensor B are 



\ Az J 



B rr B r( p B rz 

B(pr B<jxf> B(f>z 
B Z r B Z( h B zz 



U- Ay 

\ Az J 



(28) 



U • 



J-'xx D xy J-'xz 

Byx Byy By Z 

y B Z x B Z y B zz j 



(29) 



where 



U 



x/r y/r 
— y/r x/r 
1 



(30) 



is a unitary matrix. 

The NS constitutive equations yield 

P rr P(f>4> Bzz Vi 



B rz = -V 



dr ' 
dT 

q r — — k— — , 

or 



(31) 
(32) 

(33) 
(34) 
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The combination of Eqs. (25)-(27), (31)-(33) gives the following hydrody- 
namic equations: 



^ = 0, 
dr 



-id 



du. 



r — \rr] 



dr \ dr 



= -pg, 



—r\r 



'du 2 



dr \ dr I 1 \ dr 



(35) 
(36) 

(37) 



In contrast to what happens in the channel case, Eqs. (8)-(10), it is not pos- 
sible to obtain the explicit solution to the hydrodynamic equations (35)-(37), 
even with the help of an auxiliary space variable. On the other hand, the so- 
lution can be recursively found as a series expansion in powers of g. To fourth 
order the result is 



/ \ Po9 2 

u z (r) =u - - — r 

4^o 



p(r) = po, 

(4-3a)pfo 2 

i 

576r} K T 



1 + 



T ( r ) = To _ J^r 4 



Q4r] K 



1 + 



;il-9a)pgff 2 ; 
768r) K T 



+ 0(g 5 ), 
+ 0(g% 



(38) 
(39) 

(40) 



where we have taken into account that u z and T must be finite at r = 0. 
The subscript now denotes quantities evaluated at r = 0. Also, we have 
assumed that i] oc T 1_a , k oc T 1 ~ a , which corresponds to repulsive interaction 
potentials of the form [15] tp(r) oc r _/3 with a = 1/2—2/(3. The cases a = and 
ol — \ correspond to Maxwell molecules {(5 = 4) and hard spheres (f3 — > oo), 
respectively. The corresponding fluxes are 



/ \ Po9 

An = — r 



pW 



192?7o«o7o 



0(g% 



(41) 



q r (r) 



Pl9 2 r s 
16?7o 



(5 - 3a)p 2 g 2 
384^/toTo 



1 + 



+ 0(g 6 ). 



(42) 



For arbitrary g the equation of state is not as simple as in the planar case, 
Eq. (24). Elimination of r between Eqs. (39) and (40) yields 

T . T _ £L (U , _ ^ _ (l^£M (Mo _ Uzf + 0(g% (43) 



8 



3 Kinetic theory description of the channel Poiseuille flow. A sum- 
mary 



The Poiseuille flow induced by an external force in the channel geometry has 
been analyzed in the framework of kinetic theory [3,5,8-11], as well as by 
numerical simulations of the Boltzmann equation [6] and molecular dynamics 
simulations [9]. The emphasis in these papers, in contrast to that of other 
works [2,4], was put on highlighting the limitations of the NS hydrodynamic 
description (see Sec. 2.1) when the strength of the external field g is not small 
enough. 

In this Section we briefly summarize the main results derived in Ref. [5] from 
the BGK model of the Boltzmann equation. The BGK kinetic equation reads 
[17] 

% + v • V/ + m" 1 !- • (F/) = -!/(/ - / LE ), (44) 



where /(r, v,t) is the one-particle velocity distribution function, F is an ex- 
ternal force, u(t, t) is an effective collision frequency and 

f m [v - u(r,t)] 2 l 

exp {" kfb, } (45) 



/ LE (r, v,t) = n(r,t) 



m 



2Trk B T(r,t) 



is the local equilibrium distribution function. Here, m is the mass of a particle, 
ks is the Boltzmann constant, n(r, t) is the local number density, u(r, t) is the 
local flow velocity and T(r, t) is the local temperature. These hydrodynamic 
fields are defined as velocity moments of / by 

n = Jdvf, (46) 
nu = J dvvf, (47) 
B T = ^fdvV 2 f, (48) 



nk 



where in the last equation we have introduced the peculiar velocity V = v — u. 
The fluxes of momentum and energy are characterized by the pressure tensor 



P 

1 *3 



mjdvViVjf (49) 
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and the heat flux vector 



q=y/^v\/ 2 V/. (50) 



The trace of the pressure tensor is 3p, where p = nk B T is the (local) hydro- 
static pressure. The collision frequency v is proportional to the density and 
its dependence on the temperature changes in accordance with the interaction 
potential considered. For instance, v oc nT 1 ^ 2 for hard spheres, while v oc n 
for Maxwell molecules. 

For the steady Poiseuille flow in a channel, the BGK equation (44) reduces to 

( v 4 +9 t) f= -^- f ^ (5i » 



which must be complemented with the appropriate boundary conditions at 
y = ±H. On the other hand, we assume that the separation between the 
plates is large enough to allow for the existence of a bulk region —H + 5 < y < 
H — 5, where 5 is the width of the boundary layers and comprises a few mean 
free paths. Inside the bulk region the solution to Eq. (51) is expected to be 
rather insensitive to the details of the boundary conditions and depend on y 
through a functional dependence on the hydrodynamic fields. Such a solution 
was obtained in Ref. [5] (for Maxwell molecules) by means of a perturbation 
expansion in powers of g. Here we quote the hydrodynamic fields through third 
order: 



p(y) = Po 



mg \ 2 2 



+ 0(g 4 ), (52) 



/ \ Po9 2 

u z (y)=u -—y 



i + J l9 \ y 4 



30r] K, T( 



mg 
k B T 



y 2 + C 







PoVo9 2 



Po 



+ 0{g% 



T(y) = To 



l2r] K T 



y 4 + Cr{ 



mg \ 2 



+ 0(g% 



(53) 
(54) 



where ( p = § = 1.2, ( u = ±§ = 6.08, C = W = 218 - 96 > and Ct = i = 0.76. 
The BGK predictions (52)-(54) have been confirmed by an exact solution of 
the Boltzmann equation for Maxwell molecules [8] (( p — |, & — 1.0153), as 
well as by approximate solutions of the Boltzmann equation for hard spheres 
by a 13-moment method [9,11] (( p = |, ( T = H = 0.56) and a 19-moment 
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method [11] (( p ~ 1.214, (t — 0.99). Comparison with the NS predictions, 
Eqs. (14), (20) and (21), shows that the latter already fail to second order 
((NS = ^NS = ) and tQ third Qrder ^ns = ^NS = o) in g m According to Eqs. 

(52)-(54), the pressure increases parabolically from the midpoint (( p > 0) 
rather than being uniform, the velocity profile has an enhanced quadratic co- 
efficient ((' u > 0) plus a new quartic term ((„ > 0), and the temperature 
has a positive quadratic term ((t > 0). The latter is responsible for the fact 
that d 2 T/dy 2 \ y=0 > 0, in contrast to the NS prediction (12), so the tempera- 
ture presents a local minimum at y = rather than a maximum. This min- 
imum is surrounded by two maxima located at y — y max = ±^6(t^o, where 
£ Q = {jiqKqTq) 1 ! 2 /po is a reference mean free path [11]. The relative differ- 
ence between the maxima and the minimum is (T max — T )/T = 3(f.(£o/^o) 2 , 
where ho = (kBT /m)/g is the so-called scale height [18], i.e. the charac- 
teristic distance associated with the external (gravity) field. This surprising 
bimodal form of the temperature profile is an effect going beyond the Burnett 
description [10] and has been confirmed by Monte Carlo simulations of the 
Boltzmann equation for hard spheres [6,10]. The NS equation of state (24) is 
now augmented by an extra term: 



In addition to the hydrodynamic profiles, the momentum and heat fluxes are 
obtained from the kinetic theory description. In the case of the BGK model, 
the results are [5]: 




3 k q 



(u - U z f + 



(tT q 
Cp Po 



(p-p ) + O(g 4 ). 



(55) 




(56) 




(57) 




(58) 




(59) 



Q y (y) 




y 3 + o(g 4 ), 



(60) 
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29py 
12i] k, T 



y 4 +0(g 5 ) 



(61) 



As expected, the fluxes differ from the NS results, Eqs. (4), (7), (22) and 
(23). The main deviations of the hydrodynamic and flux profiles from the 
NS predictions occur for distances on the scale of the mean free path, i.e. in 
the regime where a hydrodynamic description is not expected to hold. For 
instance, the extra terms appearing in Eq. (53) are, relative to the g 2 -term 
of Eq. (20), of orders (y/£ )~ 2 and (y/io)' 4 - In addition, the ratio between 
the component of the heat flux parallel to the flow direction (q z ) and the 
component parallel to the thermal gradient (q y ) is of order (y / £o)~ 3 (ho / £o) . 
Thus, the NS description applies in the regime (y/£o) ^ (^o/^o) 1 ^ 3 ^ 1- Cm 
the other hand, the kinetic theory description, while limited here to weak 
fields, i.e. (h /£ ) 3> 1, is still valid for y ~ £ . 



4 Kinetic theory description of the pipe Poiseuille flow 

Now we are going to analyze the solution of the BGK equation for the steady 
Poiseuille flow problem in a cylindrical geometry. In that case, Eq. (44) be- 
comes 



Note that the derivative d/d(f) is understood at constant (v x ,v y ), not at con- 
stant (v r ,v<j)). In fact, dv r /d<f> = and dv^/dcf) = —v r , and so 



The exact conservation equations (25)-(27) can be reobtained from Eq. (62) 
by multiplying both sides by v r , v z and V 2 , respectively, and integrating over 
the velocity. 

As in Sec. 2, we denote by a subscript those quantities evaluated at the axis 
of the pipe (r = 0). Thus, Vq = (/c^To/m) 1 / 2 is a thermal velocity, Ao = v /u 
is a mean free path and ho = v 2 /g is a characteristic length associated with 
gravity (scale height). Since in the BGK model rj = p/u and k = 5/cs?7/2m, 
one has £ = A , where £ was introduced in the previous Section. Without 
loss of generality, we will assume a reference frame stationary with the flow 
at r = 0, so u Q = 0. We next introduce dimensionless quantities as 




(62) 




(63) 



v* = v V, r* = A 1 r, f* = n x v^f, 



(64) 
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p*=Po X p, u* = v^u, T* = T 1 T, 
u* = u 1 iy, g* = (voiso^g = \ /h . 



(65) 
(66) 



In order to simplify the notation, the asterisks will be dropped henceforth, so 
all the quantities will be understood to be expressed in reduced units, unless 
stated otherwise. 

The objective now is to find the solution to Eq. (62) as an expansion in powers 
of g: 



f = / (0) + f {1) 9 + + + 



(2) 2 , f (3)_3 



(67) 



where is the equilibrium distribution function normalized to p^ = 1, 
T(°) = 1. Similar expansions hold for the hydrodynamic fields: 



p = 1 + pWg 2 + pWg* + -.. 

u z = u {1) g + m (3) g 3 H , 

T = l + T^g 2 + T^g A + -- 



(68) 
(69) 
(70) 



where we have taken into account that, because of the symmetry of the prob- 
lem, p and T are even functions of g, while u z is an odd function. Insertion of 
Eq. (67) into Eq. (62) yields 



n-2 



(1 + A) f {n) = /g - Vf^ - £ i/»-™> (/<™> - /<&>) , (71) 



m=l 



where the operators A and X> are defined as 

^ = ^7^~ H 



V 



_d_ 

dv z 



(72) 
(73) 



The formal solution to Eq. (71) is 



/<»> = 

fc=0 



(n-l) 



n-2 
m=l 



(n— m) 



(74) 



This solution is not complete because appears implicitly on the right side 
through the dependence of /S on p( n \ and TH If the space dependence 
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of these quantities were known, Eq. (74) would give us f( n \ provided that 
the previous contributions {f( m \m < n — 1} are known. In order to get a 
closed set of equations for p^ n \ and T^, we must apply the consistency 
conditions 



/ rfv (/W - /&>) = o, 
fdvv, (/ W -/S)=0, 
/rfv, 2 (/<"> -/&>)= °" 



(75) 
(76) 
(77) 



4-1 First- order results 



In this case, = -« (1) P/ (0) , and so Eq. (74) yields 



(1) _ f(!) 



LE 



k=i 



(78) 



where we have taken into account that the operators A and V commute and 
that A k f {0) = for k > I. The conditions (75) and (77) are automatically 
satisfied. As for condition (76), it implies that 



/oo 
S(-^) fc « (1) / (0) = "I- 

fc=l 



(79) 



The simplest solution to Eq. (79) is expected to be of the form 

u {1) (r) = u 12 r 2 . (80) 



To confirm this, let us express the operator A in Cartesian coordinates, i.e. 
A = v x d/dx + v y d/dy. Consequently, only the terms with k < 2 contribute in 
Eqs. (78) and (79). Insertion of (80) into (79) then gives 



Ul2 = 



(81) 



The explicit expression for is simply 

f W = v 



1 1 1 

-i / 2 2\ 2 



/( 0). 



(82) 
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The nonzero components of the pressure tensor and the heat flux are, to first 
order, 



P W= dvv r v g fW = -r, 



1 



-1. 



(83) 
(84) 



4-2 Second- order results 



The second-order contribution to the distribution function is 

CO 

/ (2) = ftl - + E("^) fc (/S - , (85) 



fc=i 



where 



f (2) _ 



2 



//- (r: - 1) + p( 2 ) + ^ V ~ 5) 

2 



./ 



(0) 



(86) 



Equation (76) is identically satisfied, since both /f E } and P/W are even func- 
tions of v z . Conditions (75) and (77) become, respectively, 



/CO 
fc=i 



(87) 



/ dv, 2 EC-^VS = "2 / dv, z (l - A + A 2 ) 
J fc=i J 



:4+ -r 2 



(88) 



Next, by looking for a solution with a spatial dependence similar to that of 
the planar case, Eqs. (52)-(54), we write 



p (2) (r) = p 2 2r 2 , 
T^ 2 \r)=T 22 r 2 + T 2A r\ 



(89) 
(90) 



Consequently, only the terms with k < 4 contribute in Eqs. (87) and (88). 
More explicitly, 
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E(-^)*/S 

k=i 



P22 + -yCu - 5) 



-2rw r + 2(t; 2 + ^) 



^-l) + ^V- 5 ) 



- - — ( r: - ! ) + — ( r z - 0) [-4r 3 u r + 4r 2 (3^ + 7 
-24rt; r (t; r 2 + ^)+24(t; r 2 + ^) 2 ]}/( ). 

Insertion into Eqs. (87) and (88) yields 

P22 + 48T 24 = 0, 

20p 22 + 2OT22 + 2688T 24 + 12 + (80T 24 + l)r 2 = 4 + V, 



respectively. The solution is 



T, 



22 



V22 - W 



1 ToA 

50' 



160 



The explicit expression for is then 



f (2) = U-^ k fi%-vU-^ k f (1) 

k=0 k=0 



+ 



+ 



(vl - 1) 

1 



3 , 9 9\ "^2 

1 - g( U r + v 4>) + rV r ~ 1 T 



10 

1 

32 



3 + — (v 2 -h) 
10 V ; 



2 ( V r + V l) ~ 2rV r + r 2 



10 V ; 



+4r 2 (3t; 2 + ^)-4rV + r 4 ]}/ (0) . 
The second-order contributions to the pressure tensor are 



(1) 2 184 7 2 



25 10 
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Analogously, 

ql 2) = \ J dvv*v r fW - ««PW = ^r 3 . (100) 
4-3 Third- and fourth- order results 

For n — 3, Eq. (74) reduces to 

00 

/( 3) = [ ; (3) _ p/ (2) _ ^(2) _ (1Q1) 



fc=0 



where 



f (3) 



+ (« (S) - - (1) T (2) - ^ (1) \ 2 ) / (0) ] . (102) 



Up to now, we have not needed to fix the temperature dependence of the col- 
lision frequency. This implies that to second order in the field the results are 
universal, i.e. independent of the interaction potential. On the other hand, the 
results of higher order are sensitive to the potential. For the sake of concrete- 
ness, we now consider repulsive interaction potentials of the form <p(r) oc r _/3 , 
for which the collision frequency is [17] v oc 

pT -(i-a) with a = i/ 2 - 2/(3. In 
particular, the case a = corresponds to Maxwell molecules (f3 = 4), while 
the case a — \ refers to hard spheres (f3 — > 00). For this class of potentials, 
„m =p( 2 )-(l- a )T( 2 ). 

Conditions (75) and (77) with n = 3 are identically satisfied because of sym- 
metry. The structure of and suggests that has a spatial 
dependence of the form 

U W(r) = u 32 r 2 + u u r A + u 3e r 6 , (103) 



so only the terms with k < 6 contribute in Eq. (101). Insertion into Eq. (76) 
yields 



2(1367 - 206a) , 
34560m 36 + 192m 34 + 4m 32 + — >- + (1728w 36 + 16w 34 

.ZD 



149 - 36a\ 2 / 4 -3a 

f r 2 + 36m 36 + 



50 J V 160 
whose solution is 



) r 4 = 0, (104) 



4(100 -a) 89 + 9a 4 - 3a , , 

u 32 = — , „34 = — goo-, M 36 = -^5- (105) 
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Once is determined, Eq. (101) gives the explicit form of J®. From it we 
can easily get 

= J dvVrVj^ = ^ r 3 + _L r 5 ; (106) 



4(1358 - 23a) 209 - 3a 2 15 - a 



25 



+ 



50 



•r + 



160 



(107) 



Proceeding in a similar way, higher order terms can be evaluated, but the 
algebra becomes progressively more cumbersome. Here we only quote the main 
results to fourth order in g. Equation (74) gives 



/<*> = E(-^) fc [/$ - ^/ (3) - - (2) (/ (2) - /S) 



fc=0 



(108) 



By assuming that 



p (4) (r) = p 42 r 2 + p 44 r 4 + p 46 r 6 , 
T (4) (r) = T 42 r 2 + T 44 r 4 + T 46 r 6 + T 48 r 8 , 



(109) 
(110) 



it turns out that only the terms with k < 8 contribute in Eq. (108). By 
symmetry, condition (76) is identically satisfied. On the other hand, conditions 
(75) and (77) give, respectively, 



4p 42 + 192p 44 + 34560p 46 + 192T 44 + 69120T 46 + 46448640T 48 



^8(243511 - 103801a) 



625 

+ 1105920T 48 ]r 2 + 



14(95 - 56a) 



25 



+ 16p 44 + 1728p 46 + 1728T 46 



39(1 -a) 
200 



+ 36p 46 + 6912T 48 



0, 



(111) 



4(448099 - 1773928a) 



+ 



625 

131783 5995.. 

250 
1345 - 642a 



+ 20T 42 + 1728T 44 + 546048T 46 + 40388 1984T 48 



+ 80T 44 + 15552T 46 + 8736768T 48 



400 



+ 180T 46 + 62208T 48 
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+ 



11 -9a 
. 960 



+ 320T 48 



r 6 = 0, 



(112) 



where in Eq. (112) we have eliminated p 4 2-p46 in favor of T 4 2-T 48 . The solution 
to Eqs. (Ill) and (112) is 



P42 



218083 - 11035a 
1250 ' 



P44 



653 - 176a 
2000 ' 



P46 



7 — a 
4800 



, (113) 



T42 — 
Tag = 



2501129 - 38495a 



6250 
454 + 87a 

72000 ' 



T 



18 



T, 



32057 - 663a 



44 



20000 



11 -9a 
307200 ' 



(114) 



From Eq. (108) we can now evaluate the fourth-order terms in the pressure 
tensor and the heat flux. The results are 



pW=/dVT, r 2 /< 4 > 



4(5087846 - 175355a) 
3125 



106029 - 3653a , 287 - 18a 4 3-a 6 

" r T7^\ TnonTT ' ( 115 ) 



2500 6000 19200 



P (4) _ / , v 2 f (4) _ 4(5087846 - 175355a) 
F ** ~J V " 3125 



3(106029 - 3653a) 2 _ 287 - 18a 4 7(3 - a) 6 

2500 r 1200 r + 19200 r ' ^ ^ 



?W =| dvv 2 J^ - 2^V 3 ) - M « 2 (V 2 ) - T^) 

8(5087846 - 175355a) 442191 - 25799a 2 



3125 1250 

1385 -492a 4 15 - a fi 

-r 4 + — r 6 , (117) 



2000 4800 



100 -a , 97 + 9a - 5 - 3a 7 . . 

r 3 H r 5 H r 7 . (118) 

25 2400 15360 V ; 
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It can be checked that the results for p(r), T(r), u z (r), P rr (r), P^(r), P r< j,(r) 
and q r (r) we have derived are indeed consistent with the balance equations 
(25)-(27). In fact, Eq. (26) allows us to obtain P$(r) with the result 



(5) _ 3(235119 + 2780a) 3 25079 + 1097a 5 

rz ~ 12500 r + 120000 r 

71 + 9a , 23 -9a q . , 

-I ' r H r 119) 

72000 3072000 ' v ; 



5 Discussion 



When Eqs. (80), (81), (89), (90), (94), (95), (103), (105), (109), (110), (113) 
and (114) are inserted into Eqs. (68)-(70), one gets the hydrodynamic pro- 
files predicted by the BGK kinetic model through fourth order in the field. 
Comparison with the NS predictions, Eqs. (38)-(40), indicates that the latter, 
while providing the correct values of U12, «36, T 2 ^ and T 48 , do not capture the 
pressure variation (p^ 2 s = = p^f = p 46 s = 0) or the lower-degree terms of 
the velocity (w£f = u^f = 0) and temperature (Tg s = Tf 2 s = T™ s = T 4 ^ s = 0) 
profiles. 

The results derived in the previous Section strongly support the conjecture 
that the expansion in powers of g is asymptotic rather than convergent [5]. For 
instance T = l + 0.14# 2 -4.0 x 10V + • • • and T = 1-0.41# 2 -1.1 x I0 4 g 4 + - ■ ■ 
at r = 1 and r = 5, respectively. Thus, the expansion is only useful if g is small 
enough to keep the first few terms only. As we did in the planar case, Sec. 3, 
we now give the hydrodynamic profiles through third order in real units: 



p(r) = p 



mg \ 2 2 



1 + CA ^ r 



0(g 4 ), (120) 



u z (r) =u - 



Po9„ 



'0 



1 | (4-3a)p^ 2 r4 



+c«( 



576r) K T 
">!J \ ~ ,:i //'n'/dV"' 



k B T 



+ 0(g% 



T{r) = T 



1 - 



\k B T 



64r] K T 



+ 0{g% 



(121) 
(122) 



where C P = ^, ( u = (89+9a)/200, (' u = §(100-a) and Ct = ^- The structure 
of these profiles is similar to that of the planar case, Eqs. (52)-(54). In par- 
ticular, the temperature presents a local minimum at r = and a maximum 
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at r = ^32( T £ = r max , where £ = (tiqUqToY^/pq. The relative difference 
between the maximum and the minimum is (T max — T )/T = 16C|,(£o/^o) 2 > 
where /io = {kBTo/m)/g. The distance from the center of the points where the 
temperature reaches the value T max in the pipe flow (r max ~ 2.12£ ) is very 
close to that of the channel flow (y max ~ ±2.1A£ ). On the other hand, the ef- 
fect is considerably smaller in the former case [(T max — T )/T ~ 0.31(£o/ho) 2 } 
than in the latter [(T max — T )/T ~ 1.7(£o/^o) 2 ]- The equation of state is also 
similar to that of the planar case, Eq. (55), 



T = 



rp VO l 

4k 



^-(p-Po) + <% 4 ). 



(123) 



As for the momentum and heat fluxes, the results through third order are 



P rr (r) =p 

P<p<p( r ) =Po 
P zz (r) = Po 



l _ 92/Wo^ + J_ / m 9 \ 2 r 2 



25 20 \k B T 

1 _ 92 p r] 2 g 2 + _3_ ^ \ 2 ^ 2 



25 pg 20 U B T 



184 p ^ 2 7 / 
25 p3 + 1Q ^ 



mg \ 2 
r 



Po + ToU B T 

„2„2 o / ™„ x 2 



+ ^ 4 ), 
+ ^ 4 ), 

+ <% 4 ), 



192?7o«o7o 25 \k B T 



q r (r) 



r A + 0(g 4 



(124) 
(125) 
(126) 



+ 0(g 5 ), (127) 



(128) 



5*(r) = 



2mgK,Q 



1 



4(1358 - 23a)p r] 2 g 2 209 ~ 3a / m# 



5/c s 
;i5-a)p§ff 2 
64r7 /toT 



25pjj 
+ 0(g 5 ). 



50 



(129) 



As expected, they strongly differ from the NS results, Eqs. (31), (34), (41) 
and (42), with the exception of q r , in which case the error of the NS value is 
of order g 4 . 

The non-monotonic behavior of T(r) is not only an interesting effect but also 
a counterintuitive result, given that the radial component of the heat flux 
monotonically increases with the distance from the pipe axis. Consider the 
inner cylinder r < r max . Within that region the temperature increases radially 
and yet the heat flows outwards from the colder to the hotter points! The 
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solution to this paradox lies in the dramatic breakdown of the Fourier law 
(33) within the region r < r max . Following Hess and Malek Mansour [11], a 
heuristic extension of the Fourier law can be written as 

dT 

- K — = q r - e 2 V 2 g r , (130) 

where £ is a characteristic distance of the order of the mean free path. Accord- 
ing to Eq. (130), the sign of the thermal gradient results from the competition 
between q r and its Laplacian. The simple estimate V 2 qy = r~ x d (rdq r /dr) jdr ~ 
q r /r 2 shows that dT jdr > for r < £. It is easy to check that Eq. (130), with 
£ — | r max — 0.71^0, is indeed consistent with the profiles (122) and (128). If 
one characterizes the deviation from the Fourier law by means of an apparent 
thermal conductivity coefficient defined by 

dT 

Qr — ~" K app"7^~> (131) 



then one has 

o-\ -i 

,2 



^app 



K 



r, '■ - 



+ 0(g 2 ). (132) 



The above ratio vanishes at r = 0, is negative in the interval < r < r 
diverges at r = r max , and finally tends to unity from above for r r, 



max i 

max- 



The breakdown of the Newton law is characterized by an apparent shear vis- 
cosity coefficient defined by 

P'rz = ^app ■ (133) 



The ratio between this coefficients and the NS shear viscosity is 



U pi 20 \k B Tj +U{9 ' 

C' u £l d 2 T 2(2533 - 31a) t% ( du z V 4 
5CrT dr 2 175 v 2 { dr J + {9 1 1 } 

The last line of Eq. (134) clearly shows that r] app incorporates super-Burnett 
terms. Those terms can be written in equivalent alternative forms by taking 
into account that 

d 2 T 3 Vo (du z \ 2 , ( T T d 2 p , ^ 4 



7/ 



dr 2 4k \ dr J ( p p dr 2 



+ 0(g% (135) 
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Fig. 1. Hydrodynamic profiles for the case g = Q.\vqvq, as predicted by the 
Navier-Stokes description (dashed lines) and by the kinetic theory description (solid 
lines). 

As an illustration of the corrections over the NS description provided by kinetic 
theory, we compare in Figs. 1-3 the hydrodynamic and flux profiles for the 
case g = 0.1iy (kBT /m) 1 ^ 2 [which corresponds to (/io/%) 2 = 40], as predicted 
by both descriptions when only terms through third order in g are retained. 
Although higher order terms are not necessarily negligible for that particular 
value of the field, the retained terms can be expected to be enough, at least 
at a qualitative level. The curves for u z and q z correspond to hard spheres 
(a = |), but they are practically indistinguishable from those of Maxwell 
molecules (a = 0). Except for the shear stress P rz and the radial heat flux 
q r , the kinetic theory predictions dramatically differ from the NS ones. The 
hydrostatic pressure grows quadratically rather than being uniform, the flow 
velocity decreases more rapidly than expected from the NS description, and 
the temperature exhibits a non-monotonic behavior. Normal stress differences 
appear, the flux of longitudinal momentum along the longitudinal direction 




23 




Fig. 2. Same as in Fig. 1 but for the elements of the pressure tensor. 




12 3 4 



Fig. 3. Same as in Fig. 1 but for the components of the heat flux. 

being larger than the other two normal stresses [P zz > > P„)- Finally, 
the longitudinal heat flux is not only different from zero (despite the absence 
of longitudinal gradients), but it can be even larger than the radial heat flux. 
These features are qualitatively similar to those found in the rectangular geom- 
etry (cf. Sec. 3), which have been confirmed by computer simulations [6,9,10]. 



24 



6 Conclusions 



In this paper we have solved the BGK kinetic equation for the (laminar) steady 
cylindrical Poiseuille flow fed by a constant gravity field. The solution has been 
obtained as a perturbation expansion in powers of the field through fourth or- 
der and for a general class of repulsive potentials. The results exhibit a very 
weak sensitivity to the interaction potential and strongly indicate that the 
expansion is only asymptotic. A comparison with the profiles obtained from 
the Navier-Stokes (NS) constitutive equations shows that the latter widely fail 
over distances comparable to the mean free path. At a qualitative level, the 
most important limitation of the NS description is that it predicts a mono- 
tonically decreasing temperature as one moves apart from the cylinder axis. 
In contrast, the kinetic theory description shows that the temperature has a 
local minimum (T = T ) at the axis (r = 0) and reaches a maximum value 
(T = T max ) at a distance from the center (r = r max ) of the order of the mean 
free path. In the region r < r max , the radial heat flows from the colder to the 
hotter points, what dramatically illustrates the breakdown of the Fourier law. 
Furthermore, a longitudinal component of the heat flux exists in the absence 
of gradients along the longitudinal direction. Non- Newtonian effects are exem- 
plified by the non-uniformity of the hydrostatic pressure and by the presence 
of normal stress differences. 

The above effects are similar to those previously found in the case of a rect- 
angular channel. This is a non-trivial result, since both geometries are quite 
different, as can be expected from the different mathematical structure of the 
balance equations [Eqs. (l)-(3) versus Eqs. (25)-(27)] and of the kinetic equa- 
tions [Eq. (51) versus Eq. (62)]. In the rectangular geometry the relevant space 
variable y takes both positive and negative values, while the radial variable r 
is positive definite. Also, the normal stress along the gradient direction is uni- 
form in the rectangular case (P yy = const) and non- uniform in the cylindrical 
case (P rr 7^ const). At a quantitative level, on the other hand, the deviations 
of the NS profiles from the kinetic theory ones are weaker in the cylindrical 
geometry than in the rectangular geometry. For instance, the relative differ- 
ence (T max — T )/T is about 5 times smaller in the former case than in the 
latter. 

The analysis carried out in this work can be extended to the Boltzmann equa- 
tion for Maxwell molecules, as already done in the channel case [8]. The hy- 
drodynamic profiles will still be given by Eqs. (120)-(122), but with different 
numerical values for the coefficients ( p , ( u , (' u and Ct- For hard spheres, the so- 
lution can be obtained by approximate schemes, such as the moment method 
[11]. Finally, we hope that the results reported in this paper may stimulate the 
undertaking of computer simulations of the Poiseuille flow induced by gravity 
in a pipe. 
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